1 Description et objectif du cas d’étude

La transcriptomique est l’étude des ARNm qui sont transcrits à partir d’un génome. Elle permet de quantifier le niveau d’expression d’un gène transcrit en ARNm qui est ensuite traduit en protéine. La mesure de l’expression des gènes peut se faire par la méthode de la PCR quantitative (quantitative Polymerase Chain Reaction ou qPCR). On obtient ainsi une mesure relative de l’expression d’un gène par rapport à un contrôle interne (normalisation par rapport à un/plusieurs gène(s) de référence).

Les données à analyser correspondent à des niveaux d’expression de micro-ARN (miR), soit des courtes chaines d’ARN (en moyenne 22 nucléotides) qui sont des régulateurs traductionnels: en s’appariant à un ARNm complémentaire, il va le dégrader ou réprimer sa traduction en protéine.

Des données d’expression de 19 miRs mesurées par qPCR sont disponibles et l’on vous demande de les analyser à l’aide de certaines méthodes multivariées vues durant les sessions théoriques.

Vous disposez pour cela du jeu de données data_qPCR.csv qui comporte les variables “Groupe” et “Identifiant” caractérisant chaque observation, ainsi que les données d’expression de chaque miR.

2 Instructions

  1. Analysez les données log-transformées avec les méthodes suivantes:
    • 1.1. Régression logistique
    • 1.2. Régression logistique avec sélection de variables stepwise-forward
    • 1.3. Régression logistique avec Ridge
    • 1.4. Régression logistique avec Lasso
    • 1.5. PLS-DA
    • 1.6. OPLS-DA
    • 1.7. Random Forest

Avant de faire ces analyses, créez des matrices vides pour stocker efficacement les différentes sorties suivantes pour chaque classifieur (le classifieur est mis en ligne dans la matrice): * Une matrice reprenant le RMSE, RMSE-CV et la dimensionalité du modèle * Une matrice des coefficients de chaque méthode (sauf pour le Random Forest) * Une matrice des valeurs prédites par cross-validation pour chaque observation

Conseil: créez une fonction de cross-validation de type k-fold pour calculer plus facilement le RMSE-CV.

  1. Comparez des qualités prédictives des différents classifieurs :
    • 2.1. Imprimez le tableau reprenant les RMSE(-CV) et la dimensionalité des modèles.
    • 2.2. Réalisez un graphique des coefficients obtenus pour chaque classifieur et comparez les graphiques entre eux.
    • 2.3. Représentez graphiquement les valeurs prédites par cross-validation à l’aide de line plots légendés et d’un scatterplot.
    • 2.4. Dessinez des courbes ROC présentes sur un même graphe avec légende.
    • 2.5. Calculez le cut-off optimal de chaque classifieur.
    • 2.6. Enfin, présentez une table de confusion pour chaque classifieur avec le cut-off optimal et calculez le pourcentage de précision (accuracy) de ces classifieurs.

3 Importation et data management

# chemin d'accès
data.path <- "/Users/bgovaerts/Dropbox/Partage_OMICS_COURSE/Data_stat2340/Data_qPCR"
# lecture du data frame
dataqPCR <- read.csv(file.path(data.path,"data_qPCR.csv"),sep=';',header=T)

cat("dimensions \n")
## dimensions
dim(dataqPCR)
## [1] 120  40
cat("colnames \n")
## colnames
colnames(dataqPCR)
##  [1] "Groupe"         "Identifiant"    "miR_001_5p"     "miR_002_5p"    
##  [5] "miR_003_3p"     "miR_004_5p"     "miR_005_5p"     "miR_006_5p"    
##  [9] "miR_006_3p"     "miR_007_5p"     "miR_008_5p"     "miR_009_5p"    
## [13] "miR_010_3p"     "miR_011_3p"     "miR_012_5p"     "miR_013_5p"    
## [17] "let_01_5p"      "let_02_5p"      "miR_014_5p"     "miR_015"       
## [21] "miR_016"        "Log_miR_001_5p" "log_miR_002_5p" "log_miR_003_3p"
## [25] "log_miR_004_5p" "log_miR_005_5p" "log_miR_006_5p" "log_miR_006_3p"
## [29] "log_miR_007_5p" "log_miR_008_5p" "log_miR_009_5p" "log_miR_010_3p"
## [33] "log_miR_011_3p" "log_miR_012_5p" "log_miR_013_5p" "log_let_01_5p" 
## [37] "log_let_02_5p"  "log_miR_014_5p" "log_miR_015"    "log_miR_016"
# matice d'expression
X <- as.matrix(dataqPCR[,3:21]) # matrice d'expression 
rownames(X) <- dataqPCR[,2]
obsnames=rownames(X)
X <- X[,sort(colnames(X))]
miR_names <- colnames(X)
Xcat=T # indique que les variables sont des catégories bien spécifiques et pas des longueurs d'onde par exemple.  Cela a son importance dans les graphiques de coefficients et loadings
Xlog <- as.matrix(dataqPCR[,22:40])
rownames(Xlog) <- dataqPCR[,2]
Xlog <- Xlog[,sort(colnames(Xlog))]

n = dim(X)[1]
m = dim(X)[2]

classY <- rep(0, n)
classY[dataqPCR[,1]=="malade"] = 1
cat("class factor levels \n")
## class factor levels
unique(classY)
## [1] 1 0
table(classY, dataqPCR[,1])
##       
## classY malade temoin
##      0      0     57
##      1     63      0
# modification de matrices
# centrage miRs
XC <- scale(Xlog, center = T, scale = F)
XC <- data.frame(XC)

# On centre Y pour le PLS et l'O-PLS
YC <- scale(classY, center = T, scale = F)

4 Estimation des modèles de classification

Nous allons ajuster 6 types de modèles (regression logistique SW, PLS, O-PLS et logisitique ridge et logistique LASSO, random forest) afin de prédire la classe en fonction des MIR.

Pour chacun des modèles, nous allons stocker le RMSE classique (training dataset = test dataset) ainsi que le RMSE obtenu par cross-validation k-fold. Pour certains modèles, nous allons également stocker les valeurs des coefficients.

RMSE <- function(y_true, y_fit) {
  n <- length(y_true)
  sqrt((1/n)*sum((y_true-y_fit)^2))
}

##################
# kfoldCV
##################

# fonction qui fait de la k-fold CV sur une fonction de régression prédéfinie prédisant y (FUN)
kfoldCV=function(X, Y, k, FUN = ypred_FUN, seed=100, ...) {
  # Inputs
    # X, la matrice des régresseurs 
    # Y, le vecteur des réponses 
    # k, le nombre de segments et 
    # FUN : la fonction de régression qui prédit Y 
  # Outputs : 
    # ypred.cv, les prédictions par CV

  # 1. definition de k segments aleatoires
  n=dim(X)[1]
  set.seed <- seed
  indx <- sample(1:n,n)
  nseg <- k
  sn <- ceiling(n/k) # taille des segments
  
  # 2. boucle sur les sous-echantillons
  # pour calculer y-fit par cv
  yfit.cv <- numeric()
  for (i in 1:nseg) {
    if(i < nseg) {indseg <- indx[(i-1)*sn + (1 : sn)]}
    else {indseg <- indx[((i-1) * sn + 1) : n]}
    Y.train <- Y[-indseg]
    X.train <- X[-indseg,]
    X.valid <- X[indseg,]
    yfitcv <- FUN(X = X.train, Y = Y.train, XV = X.valid, ...)
    yfit.cv[indseg] <- yfitcv
  }
  
  # 3. Sortie de la fonction
  return(yfit.cv)
  }
# Définition de paramètres de modélisation
# nombre de segments en cross validation
k <- 10
# Nombre de composantes max
NCmax=10
# Liste des modèles
modnames <- c('Logistique-SW','PLSR','O-PLS','Logistique-RIDGE','Logistique-LASSO','Random Forest')
nmodels <- length(modnames) # nombre de modèles évalués
# Preparation des matrices de résultats
# Matrice des RMSE
RMSEMat <- matrix(nrow = nmodels, ncol = 3)
colnames(RMSEMat) <- c('RMSE-train','RMSE-CV','Taille du modèle')
rownames(RMSEMat) <- modnames
# Matrice des coefficients (pour tous les modèles sauf RF)
CoefficientMat <- matrix(0, ncol = m, nrow = nmodels-1)
colnames(CoefficientMat) <- colnames(XC)[1:m]
rownames(CoefficientMat) <- modnames[-nmodels]
# Matrice des prédictions
PredMat <- matrix(nrow = n, ncol = nmodels)
colnames(PredMat) <- modnames
# Tout est mis dans une liste qui sera remis à jour modèle par modèle
listr=list(CoefficientMat=CoefficientMat,RMSEMat=RMSEMat,PredMat=PredMat)
# création d'une fonction pour stocker et imprimer les résultats de chaque modèle 
printresmod=function(i,YC,yfit,yfitCV,coef,parmod,listr,pcoef=T)
{
  # i = N° du modèle
  # YC : réponses observée
  # yfit : réponses prédites
  # yfitCV : réponses prédites par CV
  # coef : vecteur des coefficients
  # parmod paramètre de dimension ou de lissage du modèle
  # pcoef : T/F pour indiquer si des coefficient sont dispo pour le modèle (pas le cas pour les RF)
par(mfrow=c(1,2))
plot(YC,yfit,pch=20,main="Y-Yfit")
plot(YC,yfitCV,pch=20,main="Y-YfitCV",ylab="yfitCV")
# calcul des RMSE 
n=length(YC)
listr$RMSEMat[i,1]=sqrt(sum((YC-yfit)^2)/n)
listr$RMSEMat[i,2] =sqrt(sum((YC-yfitCV)^2)/n)
listr$RMSEMat[i,3] = parmod 
# Récupération des coefficients et des réponses
if(pcoef==T) listr$CoefficientMat[i,]=coef
listr$PredMat[,i]=yfit
# affichage de résultats
  pander(listr$RMSEMat[i,])
  par(mfrow=c(1,1))
if(pcoef==T){  
if(Xcat==F){plot(Xval,listr$CoefficientMat[i,],type="l",ylab="Coefficients",main=paste("Coefficients ",dimnames(coefficient)[[1]][i]),xlab=Xxlab)
 abline(h=0)}
if(Xcat==T){dotchart(coef,labels=colnames(coef),
           xlab="Coefficients",main=paste("Coefficients - ",modnames[i])) } 
}
return(listr)
}

4.1 (1) Régression logistique stepwise forward

# Ajustement du modèle de regression stepwise forward
null <- glm(classY ~ 1, data=XC,family="binomial") # modele de depart
full <- glm(classY ~ ., data=XC,family="binomial") # modele avec toutes les variables
fit.sw <- step(null, scope=list(lower=null, upper=full), direction="both",trace=0) # critère d'opimisation= AIC

pander(summary.glm(fit.sw))
  Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5942 0.3373 -1.761 0.07817
log_miR_002_5p -12.62 3.313 -3.81 0.0001387
log_miR_009_5p 13.01 3.304 3.939 8.197e-05
log_miR_011_3p 2.022 0.9884 2.046 0.04076
log_miR_010_3p -1.175 0.6565 -1.79 0.07349

(Dispersion parameter for binomial family taken to be 1 )

Null deviance: 166.1 on 119 degrees of freedom
Residual deviance: 104.1 on 115 degrees of freedom
yfit.train <- fit.sw$fitted.values
namessw=names(coefficients(fit.sw))
coef <- CoefficientMat[1,]
coef[namessw[-1]]=coefficients(fit.sw)[-1]
nc <- length(coef)-1

ypred_FUN <- function(X,Y,XV,...){
  null <- glm(Y ~ 1, data = X,family = "binomial") 
  full <- glm(Y ~ ., data = X,family = "binomial") 
  fit.sw <- step(null, scope = list(lower = null, 
                                    upper = full),
                 direction="both",trace=0)
  coef <- coefficients(fit.sw)
  
  XV <- cbind(1, XV[,names(coef)[-1]])
  yfit.cv <- as.matrix(XV) %*% coef
  yfit.cv <- exp(yfit.cv)/(1+exp(yfit.cv))
  return(yfit.cv)
}

yfit.cv  <- kfoldCV(X = XC, Y = classY, k = k, FUN = ypred_FUN)

# Sauvegardes et impression de résultats
listr=printresmod(1,classY,yfit.train,yfit.cv,coef,nc,listr)

4.2 (2) Partial least squares

PLSR par la librairie pls. On estime directement les modèles par cross validation afin de trouver le bon nombre de composantes. Ensuite on reestime ce modèle là.

# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle 
ypred_FUN <- function(X, Y, XV, ...){
  fit.plsr <- plsr(Y ~ . , data = X, ...)
  coef <- coefficients(fit.plsr)
  yfit.train <- as.matrix(XV) %*% coef
  return(yfit.train)
}
## Boucle sur le nombre de composantes possibles: 1 à NCmax
RMSE.cv <- numeric()
yfit.cv <- list()
for(i in 1:NCmax) {
  res.cv <- kfoldCV(X = XC, Y = YC, k = k, 
                    FUN = ypred_FUN, 
                    ncomp = i)
  yfit.cv[[i]] <- res.cv
  RMSE.cv[i] <- RMSE(y_true = classY, y_fit = yfit.cv[[i]])
}


# Recherche du nombre de composantes à garder

plot(1:10, RMSE.cv, type = "o", main = "RMSE cross-validation")
ncomp.opt = which.min(RMSE.cv)
ncomp.opt = max(2, ncomp.opt)
abline(v = ncomp.opt, col = 2, lty = 2)

RMSE.cv <- RMSE.cv[ncomp.opt] 
yfit.cv <-  yfit.cv[[ncomp.opt]] + mean(classY)

# Ajustement du modèle avec le nombre optimal de composantes
fit.plsr <- plsr(YC ~ . , data = data.frame(XC), 
                 ncomp = ncomp.opt)
print(summary(fit.plsr))
## Data:    X dimension: 120 19 
##  Y dimension: 120 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
##     1 comps  2 comps
## X     36.76    47.46
## YC    26.76    33.89
## NULL
yfit.train <- fit.plsr$fitted.values[, , ncomp.opt] + mean(classY)
coef <- coefficients(fit.plsr)[, , ]
nc <- ncomp.opt

# Sauvegardes et impression de résultats
listr=printresmod(2,classY,yfit.train,yfit.cv,coef,nc,listr)

4.2.1 Scores

Représentation graphique des scores pour les 2 premières composantes avec la librairie MBX

par(mfrow=c(1,1))
# graphe des scores avec la librairie MBXUCL
ScatterPlot(fit.plsr$scores[,1],fit.plsr$scores[,2], createWindow=FALSE, points_labs =obsnames, main = paste("PLS score plot for PC1 and PC2"),  color=classY, pch=classY,xlab="PC1",ylab="PC2")

4.2.2 Graphe des loadings

loadings <- fit.plsr$loadings
for (i in 1:2) {
plot(loadings[,i], type="h", xaxt="n", xlab="",
     ylab = "Loading", main = paste0("Loading ", i), lwd = 3)
axis(side = 1,  at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)
}

4.3 (3) Orthogonal partial least squares

fit.opls <- opls(x = Xlog, y = classY, predI = 1, 
                 orthoI = 4, scaleC = "center", 
                 printL = TRUE, plotL = FALSE)

# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle 
ypred_FUN=function(X,Y,XV,...) {
  fit.opls <- opls(x = X, y = Y, predI = 1, 
                   scaleC = "center" , ...)
  ropls::predict(fit.opls, XV)
}

## Boucle sur le nombre de composantes possibles
RMSE.cv <- numeric()
yfit.cv <- list()

for(i in 1:NCmax) {
  rescv <- kfoldCV(Xlog, classY, k = k, FUN = ypred_FUN, 
                   orthoI = i, printL = FALSE, 
                   plotL = FALSE)
  yfit.cv[[i]] <- rescv
  RMSE.cv[i] <- RMSE(y_true = classY, y_fit = yfit.cv[[i]])
}

# Recherche du nombre de composantes à garder
plot(RMSE.cv, type = "o", main = "RMSE cross-validation")
ncomp.opt <- which.min(RMSE.cv)
ncomp.opt <- max(2, ncomp.opt)
abline(v = ncomp.opt, col = 2, lty = 2)
yfit.cv <-  yfit.cv[[ncomp.opt]] 

# Ajustement du modèle avec le nombre optimal de composantes
fit.opls <- opls(x = Xlog, y = classY, predI = 1, orthoI = ncomp.opt, 
     scaleC = "center", printL = TRUE, plotL = FALSE)
yfit.train <- ropls::predict(fit.opls, X)
coef <- fit.opls@weightMN %*% fit.opls@cMN
nc <- ncomp.opt 

# Sauvegardes et impression de résultats
listr=printresmod(3,classY,yfit.train,yfit.cv,coef,nc,listr)

4.3.1 Scores

Représentation graphique des scores OPLS

scores_p <- fit.opls@scoreMN
scores_o <- fit.opls@orthoScoreMN

par(mfrow=c(1,1))
# graphe des scores avec la librairie MBXUCL
ScatterPlot(scores_p[,1], scores_o[,1], createWindow=FALSE, points_labs =obsnames, main = paste("O-PLS score plot for PCP and PCO1"),  color=classY, pch=classY,xlab="PCPredictive",ylab="PCO1")

4.3.2 Graphe des loadings

loadings_p <- fit.opls@loadingMN
loadings_o <- fit.opls@orthoLoadingMN

plot(loadings_p, type="h", xaxt="n", xlab="", ylab = "Loading", 
     main ="Loading prédictif", lwd = 3)
axis(side = 1,  at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)

plot(loadings_o[,1], type="h", xaxt="n", xlab="", ylab = "Loading", 
     main ="Loading ortho 1", lwd = 3)
axis(side = 1,  at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)

4.4 (4) Régression logistique avec RIDGE

# Ajustement du modèle de regression pénalisé avec toutes les données
# lambda2 => RIDGE
fit.ridge <- penalized(response = classY, 
                       penalized = XC, 
                       model =  "logistic",
                       lambda1 = 0, 
                       lambda2 = 0.1, 
                       trace = FALSE)


# Optimisation du paramètre lambda2 par cross-validation
lambda2 <- optL2(response = classY, penalized = XC,   
                lambda1 = 0,  model =  "logistic", trace=0)$lambda

# Regression logistique ridge avec lambda optimisé
fit.ridge <- penalized(response = classY,penalized = XC, model =  "logistic",
                       lambda1 = 0, lambda2 = lambda2,  trace = FALSE)

yfit.train = attr(fit.ridge,'fitted')

coef <- attr(fit.ridge,'penalized') 
nc = sum(!(coef)==0)


# Calcul du  yfit.cv 
# fonction pour obtenir les y prédits
ypred_FUN = function(X, Y, XV, ...){
  fit.penalized <- penalized(response = Y, penalized = X, 
                         model =  "logistic", trace=FALSE, ...)
  
  penalizedcoeff <- attr(fit.penalized,'penalized') 
  intercept <- attr(fit.penalized,'unpenalized') 
  coef <- c(intercept, penalizedcoeff)
  
  yfit.cv=as.matrix(cbind(1,XV))%*%coef
  yfit.cv = exp(yfit.cv)/(1+exp(yfit.cv))
  return(yfit.cv)
}

yfit.cv <- kfoldCV(X, classY, k = k, FUN = ypred_FUN,  lambda1 = 0, 
                 lambda2 = lambda2)

# Sauvegardes et impression de résultats
listr=printresmod(4,classY,yfit.train,yfit.cv,coef,nc,listr)

4.5 (5) Régression logistique avec LASSO

# Ajustement du modèle de regression pénalisé avec toutes les données
# lambda1=> LASSO

fit.lasso <- penalized(response = classY, penalized = XC, model =  "logistic",
                       lambda1 = 0.1, lambda2 = 0,  trace = FALSE)


# Optimisation du paramètre lambda2 par cross-validation
lambda1 <- optL1(response = classY, penalized = XC,   
                lambda2 = 0, model =  "logistic",
                trace=0)$lambda

# Regression logistique ridge avec lambda optimisé
fit.lasso <- penalized(response = classY, penalized = XC,  model =  "logistic",
                       lambda2 = 0, lambda1 = lambda1, trace = FALSE)

yfit.train = attr(fit.lasso,'fitted')

coef <- attr(fit.lasso,'penalized') 
nc = sum(!(coef)==0)

# Calcul du  yfit.cv 
# fonction pour obtenir les y prédits (c'est la m^ême que la précédente)
ypred_FUN = function(X, Y, XV, ...){
  fit.penalized <- penalized(response = Y, penalized = X, 
                         model =  "logistic", trace=FALSE, ...)
  
  penalizedcoeff <- attr(fit.penalized,'penalized') 
  intercept <- attr(fit.penalized,'unpenalized') 
  coef <- c(intercept, penalizedcoeff)
  
  yfit.cv=as.matrix(cbind(1,XV))%*%coef
  yfit.cv = exp(yfit.cv)/(1+exp(yfit.cv))
  return(yfit.cv)
}

yfit.cv  <- kfoldCV(X, classY, k = k, FUN = ypred_FUN, 
                 lambda1 = lambda1, lambda2 = 0)


# Sauvegardes et impression de résultats
listr=printresmod(5,classY,yfit.train,yfit.cv,coef,nc,listr)

4.6 (6) Random Forest

fit.rf <- randomForest(formula=as.factor(classY) ~ ., data=Xlog,  ntree=500, mtry=3)

# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle 

ypred_FUN <- function(X,Y,XV,...) {
  fit.rf <- randomForest(formula=as.factor(Y) ~ ., data = X, ntree = 500, ...)
  yfit.cv <- predict(fit.rf, newdata = XV,  type = 'class')
  yfit.cv
}

## Boucle sur le nombre de composantes possibles
RMSE.cv <- numeric()
yfit.cv <- list()

j <- 1
mtry_tested <- seq(3,19, 2)
for(i in mtry_tested) {
  rescv <- kfoldCV(Xlog, classY, k = k,  FUN = ypred_FUN, mtry = i)
  yfit.cv[[j]] <- rescv
  yfit <- as.numeric(yfit.cv[[j]]) - 1
  RMSE.cv[j] <- RMSE(y_true = classY, y_fit = yfit)
  j <- j + 1
}



# Recherche du nombre de composantes à garder

plot(mtry_tested,RMSE.cv, type = "o", main = "RMSE cross-validation")
mtry.opt <- mtry_tested[which.min(RMSE.cv)]
abline(v = mtry.opt, col = 2, lty = 2)



yfit.cv <- as.numeric(yfit.cv[[which.min(RMSE.cv)]])-1

# Ajustement du modèle avec le nombre optimal de mtry
fit.rf <- randomForest(formula=as.factor(classY)
                       ~ ., data = Xlog, 
                       ntree=500, mtry=mtry.opt)

yfit.train <- predict(fit.rf, newdata = Xlog, type='class')
yfit.train <- as.numeric(yfit.train)-1
nc = m


# Sauvegardes et impression de résultats
listr=printresmod(6,classY,yfit.train,yfit.cv,coef,nc,listr,pcoef=F)

5 Evaluation des qualités des classifieurs

5.1 Matrice des RMSE en estimation et k-fold-validation

cat(c("Taille des segments", k))
## Taille des segments 10
pander(listr$RMSEMat)
  RMSE-train RMSE-CV Taille du modèle
Logistique-SW 0.3829 0.443 18
PLSR 0.406 0.4301 2
O-PLS 0.5736 0.4473 2
Logistique-RIDGE 0.399 0.4442 19
Logistique-LASSO 0.4174 0.458 5
Random Forest 0 0.5477 19

5.2 Plots des coefficients

for(i in 1:(nmodels-1)) {
  dotchart(listr$CoefficientMat[i,],labels=colnames(listr$CoefficientMat),
           xlab="Coefficients",main=paste("Coefficients - ",modnames[i]))
  }  

5.3 Lines plots of predictions

# Line plot to the Predictions
par(xpd = TRUE, mar = c(2,2,2,7))

plot(1:n,listr$PredMat[,1],lty=1,type="l",
     main="Predictions ordered by observation ID",xlab="Observation ID",ylab="Prediction",
     ylim=c(min(listr$PredMat),max(listr$PredMat)), col=1)  

for(i in 2:nmodels){
  lines(1:n,listr$PredMat[,i],lty=i,col=i)
}

legend("topright",1,modnames,col=1:nmodels,lty=1:nmodels,cex=0.8, 
       inset =c(-0.18,0))
abline(v=92)

Sur ce premier graphique, l’axe des x représente l’ordre dans lequel les mesures ont été effectuées. la droite verticale représente un changement dans le kit utilisé pour mesurer les expresisons des miRs. On peut constater que cela a un impact non négligeable sur les prédictions en cross-validation.

par(xpd = TRUE, mar = c(2,2,2,7))
idx <- order(classY)

plot(1:n,listr$PredMat[idx,1],lty=1,type="l",
     main="Predictions ordered by class",xlab="Observation ID",ylab="Prediction",
     ylim=c(min(listr$PredMat),max(listr$PredMat)), col=1, xaxt="n")  
axis(side=1, at = 1:n, labels = idx, las= 2, cex.axis=0.8)
for(i in 2:nmodels){
  lines(1:n,listr$PredMat[idx,i],lty=i,col=i)
}

legend("topright",1,modnames,col=1:nmodels,lty=1:nmodels,cex=0.8, inset =c(-0.18,0))
abline(v=57)

L’axe des x représente un vecteur ordonné de la classe des observations. La ligne verticale délimite les témoins des malades.

5.4 Scatter Plot matrix of predictions:

# Scatter Plot matrix of predictions
pairs(listr$PredMat,main="Scatter Plot Matrix of predictions")

5.5 Courbes ROC

# Courbes de ROC
cutoff <- seq(min(listr$PredMat)-0.1,max(listr$PredMat)+0.1,by=c(0.01))
nco <- length(cutoff)
truepositive <- matrix(nrow=nco,ncol=nmodels)
falsepositive <- matrix(nrow=nco,ncol=nmodels)
npos <- sum(classY==1)
nneg <- sum(classY==0)
listpos <- (classY==1)
listneg <- (classY==0)

for(i in 1:nco){
  for(j in 1:nmodels) {
    # Rate of true positives in the set of positive 
    truepositive[i,j] <- sum(listr$PredMat[listpos,j]>=cutoff[i])/npos
    # Rate of false positives in the set of negatives
    falsepositive[i,j] <- sum(listr$PredMat[listneg,j]>=cutoff[i])/nneg
    }
  }


plot(falsepositive[,1],truepositive[,1],type="l",main="ROC Curves Q-PCR Models",
     xlab="False Positive Rate",ylab="True Positive Rate")

for(j in 2:nmodels){
  lines(falsepositive[,j],truepositive[,j],type="l",col=j,lty=j)
  }
legend(0.6,0.75,dimnames(listr$PredMat)[[2]],col=1:nmodels,lty=1:nmodels)

5.6 Recherche du cut-off optimal pour chaque méthode

# Recherche des seuils optimaux qui minimisent la distance au point (0,1)
cutoffopt=listr$PredMat[1,]
for(i in 1:nmodels){ 
  distances2 <- falsepositive[,i]^2+(1-truepositive[,i])^2
  cutoffopt[i] <- cutoff[which(distances2==min(distances2))][1]
  }  
pander(t(cutoffopt))
Table continues below
Logistique-SW PLSR O-PLS Logistique-RIDGE Logistique-LASSO
0.5572 0.5472 0.6272 0.5472 0.5472
Random Forest
0.007205
# Dessin des cutoff sur un scatter plot
for(i in 1:nmodels){ 
  plot(YC,listr$PredMat[,i],main=paste("Cut-offs - Method:",
                 dimnames(listr$PredMat)[[2]][i]),ylab="Prediction")
  
  abline(h=cutoffopt[i])
  } 

5.7 Matrices de confusion avec les cut-offs optimaux

PredictClass <- matrix(nrow = n, ncol = nmodels)
dimnames(PredictClass) <- dimnames(listr$PredMat)
TrueClass <- rep(1,n)

for(i in 1:nmodels){
PredictClass[,i] <- listr$PredMat[,i]>=cutoffopt[i]
}

# Calcul les nombres de bien classés avec les seuils optimaux
ClassTable <- matrix("",nrow = nmodels * 3, ncol = 2)
dimnames(ClassTable)[[1]] <- rep(c("","Class0","Class1"),nmodels)

dimnames(ClassTable)[[2]] <- c("FALSE","TRUE")
dimnames(ClassTable)[[1]][((1:nmodels)*3)-2] <- modnames

for(i in 1:nmodels){
  ClassTable[(i-1)*3+(2:3),] <- table(classY,PredictClass[,i]) 
}

cat("Classification Tables")
## Classification Tables
pander(ClassTable)
  FALSE TRUE
Logistique-SW
Class0 42 15
Class1 12 51
PLSR
Class0 39 18
Class1 11 52
O-PLS
Class0 38 19
Class1 10 53
Logistique-RIDGE
Class0 43 14
Class1 13 50
Logistique-LASSO
Class0 42 15
Class1 16 47
Random Forest
Class0 57 0
Class1 0 63
# accuracy 
pc.accuracy <- c()
for(i in 1:nmodels){
  pc.accuracy[i] <- sum(diag(table(classY,PredictClass[,i])))/n
  }

pc.accuracy <- as.matrix(pc.accuracy)
colnames(pc.accuracy) <- "pc.accuracy"
rownames(pc.accuracy) <- modnames
pander(pc.accuracy)
  pc.accuracy
Logistique-SW 0.775
PLSR 0.7583
O-PLS 0.7583
Logistique-RIDGE 0.775
Logistique-LASSO 0.7417
Random Forest 1